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We study a gas of repulsively interacting bosons in the first excited band of an optical lattice. 
We explore this p-band physics both within the framework of a standard mean-field theory as well 
as with the more accurate generalized Gutzwiller ansatz. We find the phase diagrams for two- and 
three-dimensional systems and characterize the first Mott-states in detail. Furthermore, we find that 
even though the p-band model has strongly anisotropic kinetic energies and inter-flavor interaction 
terms are missing in the lowest band theory, the mean-field theory becomes useful quite rapidly 
once the transition from the Mott-insulator to the superfluid is crossed. 
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I. INTRODUCTION 

Systems of cold atoms in optical lattices have seen a 
dramatic experimental progress in the recent past [l], Q • 
Due to realization of optical lattices, low densities, and 
low temperatures, a fantastic degree of control has been 
obtained which has made detailed studies of strongly 
correlated quantum systems possible. For example, the 
Mott-superfluid transition [H, 0] has been successfully ob- 
served in optical lattices. This transition, due to the 
competition between kinetic energy and repulsive on-site 
interactions between lowest band bosons, can occur even 
at T — and is therefore driven by quantum fluctua- 
tions. For large interactions, the energy is minimized in 
an incompressible state with fixed atom numbers at each 
lattice site, while for weaker interactions the kinetic en- 
ergy favors atomic tunneling which drives the system into 
a superfluid. 

The early experiments were confined to the lowest 
band and while increasing interactions can make excited 
band populations non-negligible Q , the lowest band still 
dominates. In fact, for very strong interactions, it has 
been theoretically shown that the lowest band Mott in- 
sulator turns into a Mott insulator at the p-band Q. 
Experimentally, however, atomic population residing on 
the excited bands is obtained by couple atoms from the 
lowest band to the excited bands. This was experimen- 
tally demonstrated by accelerating the lattice for a short 
period or more recently by coupling atoms from the 
lowest band Mott insulator into the first excited p-band 
of the lattice via Raman transitions between bands [g] . In 
the latter of these two, it was in particular found that the 
lifetimes of p-band atoms are considerably longer than 
the tunneling time-scale in the lattice and they were also 
able to explore how coherence on the excited band es- 
tablishes. These experiments pave the way to explore 
also equilibrium physics of the purely p-band bosons @ 
and furthermore provide possible routes to realize super- 
solids [k| or novel phases [HI E3] on the excited bands 
of an optical lattice. An alternative way to populate 
higher bands is by considering fermions with a filling fac- 
tor larger than one [H, 0, In this case the Pauli 
exclusion principle ensures that the fermions that can- 



not populate the lowest band, must occupy the excited 
bands [lq |. 

In this paper we explore the properties of bosons oc- 
cupying the first excited bands of an optical lattice. In a 
periodic potential where the lattice depths are equal in 
all directions, the (non- interacting) bands are degener- 
ate and a multi-flavor description of the quantum states 
of atoms is required d, [ijj • This fact together with the 
non-isotropic tunneling on the p-band add new features 
and possibilities both for the description of the superfluid 
as well as insulating phases. For example, onsite flavor 
changing collisions can induce fluctuations in the num- 
ber of atoms of different flavors even in the insulating 
phases, giving them non-trivial characteristics. Further- 
more, such collisions together with anisotropic tunneling 
cause different types of phaselockings (both locally as 
well as between sites) between flavor condensates in the 
broken symmetry phases. 

While in many places we confirm the general picture 
provided by the somewhat simplified model considered 
by Isacsson and Girvin [9(. Nonetheless, we also find 
differences which arise due to; the use of real Wannicr 
functions (as opposed to the approximated ones given 
by a harmonic ansatz), through the inclusion of nearest 
neighbor tunneling in all directions, or through difference 
in accounting for the inter-flavor interactions. For fu- 
ture reference, we also compare the Gross-Pitaevskii type 
mean-field theory with the more accurate Gutzwiller ap- 
proach and find the parameter regions where the Gross- 
Pitaevskii approach is reasonably accurate. 

The paper is organized as follows. In Sec. |TT] we de- 
rive our model Hamiltonian and by taking anharmonic- 
ity of the lattice potential into account, we outline un- 
der what conditions the physical description can be re- 
stricted to the first excited p-band. We then proceed 
by deriving mean-field Gross-Pitaevskii equations for the 
p-band bosons and discuss salient features of their solu- 
tions for a homogeneous system both for two-dimensional 
as well as for three-dimensional systems in Sec. IIIII In 
Sees. II VI and IIVBI we study the p-band physics in two- 
and three-dimensional systems employing the Gutzwiller 
ansatz and outline the ways how these solutions differ 
from the mean-field ones. We conclude with a brief dis- 
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II. FORMALISM 

The microscopic Hamiltonian for the dilute Bose gas 
at low temperatures in a trap is given by 



H r , 



drip' (r) 
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(1) 
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where fi is the chemical potential, m the atomic mass, 
g is the interatomic interaction strength, and ip{r) and 
ip' (r) are the bosonic annihilation and creation operators 
respectively, while V(r) is the external trapping potential 
which in this work is taken to be a lattice potential 



V(r) = V L 
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where d is the lattice spacing and Vl the lattice depth. 
For a deep lattice it is reasonable to expand the field 
operators in terms of the localized Wannier functions. 
Here we go beyond the usual lowest band Hubbard model 
by also including the first excited states (p-band). In a 
three dimensional lattice this implies an expansion of the 
field operators 



■0(r) = ^2w aA (r)ip 



(3) 



where i = {ix,iy,iz) labels the lattice site and a G 
{0, x, y, z} is the flavor index. The bosonic operators ty a ,\ 
annihilate a boson of flavor a from the site i. We compute 
the Wannier functions from the ideal gas band structure 
calculations. In this paper we assume that the system has 
been prepared on an excited p-bands and in the following 
set the population of the lowest band to zero. 

Substituting the operator expansions into the Eq. ([1]) 
and ignoring all but the leading order onsite interactions 
and nearest neighbor tunneling processes we derive our 
fundamental Hamiltonian 

H = H + H nn + H FD , (4) 

where the ideal part is given by 

Hq = Y^ -W^iVv.i ~ E E ^V^iVvj- (5) 

i a, a <i,j>„ 

Here X)<ij> indicates the sum over nearest neighbors 
in the direction a € {x, y, z}. Since the Bloch functions 
diagonalize the single particle Hamiltonian, there are no 
interband hopping terms in the Wannier representation 



considered here The terms originating from inter- 

atomic interactions are given by 



Hnn = E E "^P^ 1 ~~ 1 ) 
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and 



i era' ,<j^<j' 



(7) 



+ $l> i'/'l'.i^i^i), 



where Hp d contains terms that describe flavor chang- 
ing collisions which transfer atoms between bands. This 
term has a formal similarity with terms responsible for 
spin-dynamics in a spinor condensates []~9L [20j| - However, 
the strength of these terms is comparable to other in- 
teraction terms as opposed to spinor condensates where 
it is usually small, being proportional to the difference 
between singlet and triplet scattering lengths (for spin-1 
spinor condensate). 

It should be kept in mind that there are circumstances 
when nearest neighbor interactions [To| or particle as- 
sisted tunneling processes [21( might give rise to new 
physics. These contributions are not included in the for- 
mulation presented here where our focus is in the most 
typical parameter regimes. 

The various coupling strengths in the lattice model are 
related to g through 



g / dvw a:i (v) 2 w a - :i {vf 



and the tunneling coefficients are given by 



(8) 



drw a i(r) 
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V(r) 



w CT 4 + i Q (r), (9) 



where by i + l a we indicate the neighboring site of i 
in the direction a. When the lattice is symmetric, 
the tunneling strength on the lowest band is indepen- 
dent of direction. However, this is not true for the p- 
band where the directional dependence of the tunneling 
strength must be kept, as the overlap integrals are very 
different depending on whether one is integrating along 
the node of the Wannier function or orthogonal to it. 
This indeed has im por tant consequences for the physics 
in these systems (9l. 1 1 lL HEj 

It should be further noted, that since the parameters 
of our model are computed using real Wannier functions, 
we find some, not only quantitative, but also qualitative 
differences from the commonly used models build on the 
harmonic approximation. In particular, many degenera- 
cies appearing in the harmonic approximation are absent 
when real parameters are used. 
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A. Validity of p-band single-band approximation 

In a harmonic potential, two atoms on the first excited 
states have an energy 2 X 7ia;(3/2 + 1). This is equal to 
the energy of one atom on the ground state and one atom 
on the second excited state. This suggest that collisions 
between p-band atoms can populate also the lowest s- 
band and the (i-bands. This would clearly restrict the 
validity of the models residing purely on the p-bands. 

However, a real site in an optical lattice is not ex- 
actly harmonic and this anharmonicity implies that the 
above processes are normally off- resonant. The devia- 
tion between the real lattice potential and the harmonic 
approximation is given by 



AV 



V L [sin 2 (nx/d) + sin 2 (ny/d) + sin 2 (wz/d)] 



-Vltt 2 



( 10 ) 

In the first order perturbation theory around the har- 
monic approximation, we find that in the limit of deep 
lattices (Vl/Er ^> 1, where Er is the recoil energy) the 
detuning 2£'i i o,o — ^o,o,o — ^i,i,o ; where subsripts de- 
note quantum numbers u x , u y , and v z of the harmonic 
oscillator states, vanishes. This would be related to a 
process where two atoms from the p-band scatter into 
one ground state atom and one atom occupying the state 
\v x — l,i/y — l,v x — 0). Even though this detuning re- 
mains zero at first order in anharmonicity, the process 
has a vanishing matrix element and can therefore be ig- 
nored. 

On the other hand, a process where two atoms from the 
p-band scatter into one ground state atom and one atom 
on the state \v x — 2 1 v y = 0, v z = 0) (for example) can oc- 
cur. For this process the detuning 2£ , 10 ,o — -E'o,o,o — -^2,0,0 
approaches a constant value of — 2/3 Er in the limit of 
deep lattices. Note that the oscillator energy Ku> has a 
\fVh dependence in the same limit, so even though the 
detuning approaches a constant for deep lattices, it be- 
comes small relative to the harmonic oscillator energy 
scale. From this we can conclude that as long as the 
bandwidths and interactions are very small compared 
to recoil energy, we can safely ignore d-band atoms and 
processes which would scatter atoms from the p-band to 
other bands. 

We also note that one way to prevent atoms on the_ 
band to populate the s-band was outlined in Ref. 
Here, fermionic atoms are occupying the lowest band 
and due to atom- atom interactions, the p-band atoms 
are blocked from occupying the lowest band. 



III. GROSS-PITAEVSKII APPROACH 

In the mean-field approach we replace the operators 
ipa.i with complex numbers ip a ,i- This approximation 
amounts to a coherent state ansatz in each site. In a 



Fock representation this is given by 



(11) 
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where ipi >a — (V'i.a) is the order parameter for the flavor 
a at site i = (i x ,i y ,i z )- This mean-field approximation 
is expected to be reasonably accurate in the superfluid 
phase when interactions are much weaker than the tun- 
neling strengths. In this same regime the effects due to 
the (i-band atoms, can also be safely ignored as long as 
the tunneling strengths are much smaller than the anhar- 
monicity induced detuning discussed earlier. Using the 
coherent state ansatz we can derive the equations of mo- 
tion for the order parameters ip a from the Euler-Lagrange 
equation 



dL 



d dL 



fyta dt {d^J 

with the Lagrangian given by 



0. 



(12) 



Hmf- 



(13) 



Here Hm f is the mean-field approximation for the Hamil- 
tonian in terms of the coherent state amplitudes. What 
we find are the discretized versions of the Gross-Pitevskii 
equation for each flavor. These equations are non-linear 
and coupled, but can be solved numerically without too 
much difficulty. Furthermore, in some special cases ana- 
lytical results can even be derived. We choose the lowest 
band tunneling strength as our unit of energy and lat- 
tice spacing as our unit of length. Then, for a three- 
dimensional lattice, the Gross-Pitaevskii equations for 
different p-band flavors read 



ih 



dt 



ih 



dt 



- J2a tx > a bPi+l a ,x - 2-0i,x + 1pi-l a ,x] 

+ [9xx\lpi,x\ 2 + ^9xy\^i, y? + ^9xz\A,z\ 2 ] 4>i,x 

(14) 

+ [gyy\lPly\ 2 + ^9xy\^i,x\ 2 + 2g V z\lpi,z\ 2 } i>i,y 

+^ixn y + ^iMy> 

(15) 
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and 



ih 



at 



+ [9zz\lpiA 2 + 2 9xz\lpux\ 2 + 2 9yz\i>i, y \ 2 ] -0i,z 

(16) 

In these equations the first term on the right hand side is 
due to the kinetic energy in the lattice, the second term 
originates from the density-density interactions, while 
the last terms are due to the flavor changing collisions. 
The generalization for the two-dimensional system with 
only two flavors is straight forward. 



A. 2-dimensional lattice 

In a two-dimensional system we have two degenerate 
p-bands. On a mean-field level it is easy to investigate 
the lowest energy wavefunctions in the broken symmetry 
phase. When the lattice is very deep, the energy mini- 
mization can be done in each site separately by ignoring 
the tunneling term entirely. In this way we find that the 



lowest energy state in each site is given by ip x 



a i<t> 



A/2 



and i/jy 
wavefunction 



i<f>±iT /2 j This corresponds to an onsite 



(V>(x)) = w x (x)ip x + w y (-x)ipy. 



(17) 



Since the Wannier functions are related to each other and 
can be expressed as w x (x) = f(x)u>o(x.) and w y {x.) — 
f(y)w (x), this implies 



(V>«) 



o»0 



V2 



u/ (x) (f(x)±if(y)). 



(18) 



For deep lattices the Wannier functions approach the har- 
monic oscillator states and f(x) ~ x. We can then clearly 
see that the mean-field state corresponds to a vortex or 
anti- vortex state with an angular momentum ±1 along 
the z-axis. 

Within this approximation, any configuration of either 
vortex or anti-vortex states at each site are degenerate. 
However, when the tunneling term is non-zero, phases 
of the order parameters in different sites must be corre- 
lated properly if the energy is to be minimized. When 
tunneling strengths are positive, the lowest energy con- 
densed state has the same phase at each site. However, 
on the p-band the tunneling strength for a flavor is nega- 
tive in the direction of the node in its localized Wannier 
function. In this case, the lowest energy state has a it 
phase-difference between neighboring sites. For a two- 
dimensional system it is possible to find the mean-field 
state which minimizes the onsite problem as well as the 
tunneling problem simultaneously and this state amounts 
to a checkerboard (or anti- ferromagnetic) ordering of vor- 
tices and anti-vortices. 



This is easy to see, since if at some site we have a vortex 
state ~ (x+iy) and we aim to minimize the kinetic energy 
along y-direction, then the neighboring site should have 
a same phase for the x-flavor while having a 7r-phaseshift 
for the y-flavour. This implies an anti- vortex state ~ (x— 
iy). If we then try to minimize the kinetic energy along 
x-direction, we see that the x-flavor should experience a 
7r-phaseshift, while for the y-flavour the phaseshift should 
vanish. This implies an anti- vortex state ~ e l7T {x — iy) 
with an additional overall phaseshift of it. 



B. 3-dimensional lattice 

In a three-dimensional lattice we have three degenerate 
bands, which opens up for novel phenomena not present 
in the two-dimensional case. Mimimizing the onsite prob- 
lem we find that the lowest energy configuration becomes 




1 



^e** ( exp(2«/3) 
exp(47ri/3) 



(19) 



where n<r is the total onsite atom number and (j) is a ran- 
dom phase. The onsite wavefunction with equal number 
of atoms in each flavor has an unit angular momentum 
per atom which points not along the main axes, but di- 
agonally L oc (±1,±1,±1). Again minimization of the 
kinetic energy necessitates a special ordering of angular 
momentum in each site. In the three-dimensional lattice 
the nearest neighbor angular momenta (in the direction 
a) are related by a relation L(i + e Q ) = i? ,(7r)L(i), where 
R a (ir) is a rotation of it around the axis a. 

The above results depend crucially on the magnitude 
of the inter-flavor coupling strengths g xy = g xz = g yz rel- 
ative to the magnitudes of the g aa terms. In particular, 
it only holds when g xy < g xx /3. If one approximates the 
Wannier functions with the harmonic oscillator states one 
finds that g xy = g xx /3, but when real Wannier functions 
of an ideal Bose gas are used, g xy < g xx /2> for fairly deep 
lattices and the above result holds. That said, the result 
may be different for a shallow lattice. Furthermore, it is 
also unclear what is the effect of the interaction-induced 
dressing of the Wannier functions [22j on the magnitudes 
of the effective p-band coupling strengths. If it turns out 
that under some circumstances inter-flavor coupling is 
larger and g xy > g xx /?>, then the lowest energy configu- 
ration breaks the permutational symmetry and is given 
by the vortex states 



1 

exp(±7ri/2) 




(20) 



where the angular momentum points along the z-axis. 
The vortex-anti-vortex states with angular momentum 
along other axes are degenerate with the one shown here 
explicitly. It is seen that the state (|19p has a mutual 
120° phase difference between the flavors, reminiscent of 
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three interacting spin 1/2-particles placed on the corners 
of a triangle. The state (|20[) . on the other hand, shows a 
mutual 90° phase pattern. In particular, the interaction 
terms proportional to g a /3, with a — (3, favors a 90° pat- 
tern, while those with a ^ (3 favor a 120° configuration. 

In order to achieve a better understanding how the 
particular limiting case g xy = g xx /3 comes about, let us 
again write the Hamiltonian as H — H nn + Hpjj where 
in the mean-field approximation we have 

H nn = 9xx [n 2 x + n 2 y + n 2 ] 

+4g xy (n x n y + n x n z + n y n z ), 

(21) 

Hfd = 2g xy [cos(A xy )n x n y + cos(A xz )n x n z 
I 

+ cos((A a;z - A xy ))n y n z ] . 
M = 



Thus, we have rewritten the single site problem in a 
quadratic form for the n a variables. For the general 
case, the eigenvalues are not analytically solvable. How- 
ever, assuming g xy < g xx /i we may use the fact that we 
know that the energy is minimized for A xy = 2tt/3 and 
A xz = 4-7r/3 and we then obtain 

^i = gxx — 3g xy , 

A 2 = g xx ~ 3g xy , (24) 

A3 = 9xx + Ggxy 

Since g xx > 3g xyi the matrix M is positive definite. How- 
ever, putting g xx < 3g xy into Eq. (|24|) results in a non 
positive definite matrix and we can conclude that the 
120° phase symmetry is broken in such a case. The pos- 
sibility of the broken permutational symmetry was also 
noted in Ref. @. 

We should point out that all the results rely on having 
an isotropic lattice configuration. Any deviation from the 
symmetric lattice will break this degeneracy and give a 
preferred direction for the axis of angular momentum. 
When kinetic energy is included the ordering of vortex 
anti- vortex states between sites is the same as in the two- 
dimensional system. 

IV. QUANTUM STATES 

The Gross-Pitaevskii mean-field approximation usu- 
ally provides a sufficient description when one consid- 
ers the superfluid phase with a large number of atoms 
per site. However, the local coherent state ansatz is not 
necessarily all that good when the average onsite occu- 
pation number is small. Furthermore, the mean-field de- 
scription fails completely when the system is in a Mott 



Here, n a = \ijj a \ 2 , and A a p = 4> a — tfip, with <f> a being 
the phase of (ip a )- The energy functional can be written 
in the form 



E[(i),W,(^)]=n T Mn, (22) 



where n = (n x ,n y ,n z ) and 



(23) 

I 

insulator phase i.e. when the onsite atom distribution 
is greatly sub-Poissonian. Therefore, to accurately de- 
scribe the system properties in this regime a more pre- 
cise many-body wave function is needed. This will also 
provide insight into which parameter regimes where the 
mean-field picture is a good approximation. 

We assume that the state vector has the generalized 
form of the Gutzwiller approximation [23|. This is a 
product of on-site quantum states expanded in terms of 
the Fock states |n) of the multiple flavor system 

lV')=IIE/n ) |n)i (25) 

i n 

where the index i runs over all lattice sites. The expan- 
sion coefficient /n is the Gutzwiller amplitude of the 
particular on-site Fock state. For our purposes, in the p- 
band the relevant subspace is covered by the Fock states 
of the form |n) = \n x , n yi n z ), where for example, n x is 
the occupation number of the p K -flavor. 

Within the Gutzwiller approximation the energy of the 
system becomes a functional of the amplitudes /n ■ By 
utilizing a conjugate gradient method we minimize this 
functional giving the system ground state at zero temper- 
ature, T = 0. Several aspects regarding the minimization 
were discussed in Ref. 24]. Here we only mention that 
the sum over n must be cut off and in our numerical 
scheme we include all the states with n CT < 8 in 2D 
and n a < 6 in 3D. Throughout this section we will 
choose the lattice amplitude to be fixed and instead as- 
sume that the ratio between tunneling and onsite inter- 
action can be controlled via Feshbach resonances. The 
Wannier functions are calculated for a relatively deep 
lattice, Vl = 15-Er 



yxx 

2g xy (2 + cos(A a , !/ )) 

2g xy (2 + cos(A 2 , !y )) 

H.r.i 

2g xy (2 + cos(A l:2 )) 2g xy (2 + cos((A^ - A xy ))) 



2g xy (2 + cos(A xz )) 
2g xy (2 + cos((A X2 - A xy ))) 
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A. 2-dimensional lattice 

For the two-dimensional lattice, the onsite Fock states 
then consist of p x and p y terms only. To get insight 
about possible correlations between neighboring lattice 
sites our effective computational subspace contains four 
lattice sites with two sites in both spatial directions. The 
computational 4i/Uoo-fJ,/Uoo parameter region is chosen 
to be such that the total number of atoms per site is 
relatively small. To investigate the effects of quantum 
fluctuations, this is the region of most experimental in- 
terest and it is also favorable numerically with reasonable 
cut-offs. 

Some resulting properties can be seen in Fig. [TJ The 
absolute value of the two condensate order parameters 
KV'x)!, 1(^)1 are plotted in Fig. Q] (a). As in the stan- 
dard s-band Bose-Hubbard model, the phase-space con- 
sists of Mott insulating lobes and superfluid regions. This 
is further evidenced in (b) where the total atom num- 
ber rix is shown; within the Mott lobes nx attains an 
integer value. Due to symmetry reasons, it is not sur- 
prising that the absolute values of the two flavor order 
parameters are identical. However, in the SF phase, 
there is a phase difference of ±7r/2 between the two 
flavors e.g., (ip x ) is real while (ip y ) is imaginary. This 
suggests that the on-site ground state is a vortex; a 
result in agreement with our Gross-Pitaevskii calcula- 
tions. Indeed, a plot of the scaled angular momentum 

\L\/n T = L z .i = -i N>l ;i i> y ,i - i>l ti $x,i) /n T given in 

Fig. [1] (d) verifies that for strong tunnclings and large on- 
site atom numbers the angular momentum is quantized. 
The existence of the vortex solution is also supported by 
the work of Watanabe and Pethick [25|. Namely, in a 
single harmonic trap within the the mean field approxi- 
mation the energy functional is of the form 

Emfho = - ~ [1 - (An)] sin 2 (26) 

where 7 is the effective coupling constant, An — n x — n y 
population difference between the two flavors, and (f> is 
their relative phase. Eq. (|26j) is clearly minimized when 
An = and 4> — ±7r/2. Physically this means that 
the repulsive interaction favors a vortex solution above a 
non-vortex one. 

As discussed in the previous section, in the mean-field 
limit the vortices on the lattice tend to order themselves 
in a form of a checkerboard pattern with neighboring 
vortices and anti- vortices. According to our Gutzwiller 
results this is true also more generally in the superfluid 
phase. In fact, our ground state of anti-ferromagnetic 
like vortex ordering is similar to the staggered- vortex su- 
perfluid state discussed in Ref. [26| for a square optical 
lattice in an effective staggered magnetic field. However, 
in the p-band such a state appears even in the absence 
of effective magnetic fields. 

The physics appearing for the multi-flavor Mott insu- 
lating states is possibly even more interesting. For exam- 
ple, as seen in Fig.[T](c), onsite number fluctuations An 2 , 




FIG. 1: (Color online) Properties of the two-dimensional two- 
flavor Bose-Hubbard model as a function of the chemical po- 
tential and the inverse interaction strength 4t/Uoo where the 
factor 4 derives from the number of nearest neighbors. For 
concreteness the parameters were computed for a lattice depth 
of Vl = 15Er. The various plots show: order parameters (a), 
total atom number (b), atom number fluctuations (c), and 
angular momentum per particle (d). 

(or equivalently An 2 ) for the individual flavors are not 
necessarily zero. For the Gutzwiller ansatz wave func- 
tion (|2"5)) . no correlation between sites is allowed. As an 
outcome, for odd total number of atoms n-r there is a set 
of degenerate Mott states, e.g. with ut = 1 all onsite 
interaction terms vanish and the state \n x — l,n y = 0) 
is degenerate with \n x = Q,n y = 1) or any linear com- 
bination of these. However, tunneling between sites will 
normally break these degeneracies. The Gutzwiller ap- 
proach is not able to capture such effects and therefore 
the kinetic energy term 

f = J2 E MlAd ( 27 ) 

c,a <i,j> Q 

is taken into account within second order perturbation 
theory. We focus on the lowest Mott, m = 1, and it turns 
out that the degeneracy is indeed lifted and the ground 
state shows a anti-ferromagnetic vortex structure. We 
note that the favorability of a vortex state in the tit = 1 
Mott state relies to the non zero value of the transverse 
tunneling rate. If this tunneling is compeletely neglected 
the energy is minimized by a ferromagnetic state [9(. In 
the second lowest Mott, nr = 2, the picture is simpler 
because the interactions break the degeneracy and no 
perturbation theory is needed. The state \n x — l,n y — 
1) is favored over \n x — 2,n y — 0) and \n x — 0, n y = 2) 
due to the vanishing of the self terms proportional to 
n a ,i (?V,i - 1) in Eq. ©. 

We further illustrate our results by plotting the abso- 
lute values of the Gutzwiller amplitudes f a in Fig. [2] 
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FIG. 2: (Color online) Absolute values of the Gutzwiller am- 
plitudes at a single lattice cite. The left figure (a), shows 
the atomic distribution for a superfluid ground state with 
At/Uoo = 0.04 and fi/U a = 0.7. On the right figure (b), 
a Mott insulator state is plotted for At/Uoo = 0.01 and 
n/Uoo = 0.7. Expectedly, in this insulator phase only the 
Fock state \n x ,n y ) — jl, 1} is populated within the Gutzwiller 
approach. 



as bar graphs. Hence, these are the probabilities of the 
onsite state to be at a given Fock state \n X) n y ). In Fig. 
[2] (a) the single-site amplitudes of a superfluid state are 
given when At/U 00 = 0.04 and [i/Uqo = 0.7. This state 
is clearly a superposition of many Fock states whereas 
the Mott state of Fig. [2] (b) contains only one state. 
This MI state is the minimum energy configuration for 
At/Uoo — 0.01 and (<l/Uqo = 0.7. Due to the small num- 
ber of atoms, the superfluid atomic distribution depicted 
in Fig. [5] (a) is still sub-Poissonian. It is also evident from 
the figures that the states with extensive populations are 
negligible justifying our numerical cut-off at n a < 8. 



B. 3-dimensional lattices 

In a symmetric three dimensional lattice the p-band 
is described in terms of 3-flavors. In the Mott insulator 
with only one atom per site, for the same reason as for 
the two-dimensional case, the ground state is strongly de- 
generate within the Gutzwiller ansatz. As argued above, 
such states are not true eigenstates of our Bose-Hubbard 
Hamiltonian, and again for relatively deep lattices the 
breaking of this degeneracy, and hence the permutational 
symmetry breaking, is well described within second or- 
der perturbation theory. Using real Wannier functions 
to compute the model parameters we find that, in a the- 
ory which takes the kinetic energy into account perturba- 
tively, the ferromagnetic state where only one flavor is oc- 
cupied has a lower energy than either anti-ferromagnetic 
states with checkerboard ordering or striped phases. 

With only two atoms per site the condensate order 
parameters naturally vanish when entering the Mott in- 




FIG. 3: (Color online) The condensate order parameters for 
the three-dimensional lattice. 

sulating regime, but the local angular momentum (L) = 
-j= (±1,±1,±1) is non-zero and (L 2 ) is equal to 6. An- 
gular momentum per particle \J^2 a (L a ) 2 /nT is 1/2 in 
this state and is in a marked contrast to the superfluid 
regime, where the onsite angular momentum per particle 
is equal to one. In a superfluid phase, the half-quantum 
vortex can occur in multi-component systems and can 
be pictured as a vortex in one of the component with 
the vortex free component filling the vortex core [13, HH . 
However, being non-zero the expectation value of the an- 
gular momentum is in qualitative agreement with the 
Gross-Pitaevskii solution even in the Mott lobe. More 
explicitly, the minimum energy state in each site is max- 
imally entangled angular momentum eigenstate given by 

|V) i [e #1 |H0) + e^ 2 |101) + e^ 3 |011)] , (28) 
V3 

where the amplitudes have 27r/3 phase-differences. 

For three atoms per site, the lowest energy Mott insu- 
lator state has the onsite wavefunction ip = | 111) , which 
was also found for the corresponding state in the two- 
dimensional lattice. Importantly it should be noted, that 
commonly used harmonic approximation for the Wannier 
states predicts the properties of (for example) this insu- 
lating phase incorrectly If harmonic oscillator states are 
used to approximate Wannier wavefunctions, the insulat- 
ing state with 3 atoms per site is degenerate with more 
complicated superposition states, but these degeneracies 
are removed once real Wannier states are used to eval- 
uate the parameters of the theory. As we have demon- 
strated, tunneling between sites will remove the onsite 
degeneracies among the Mott insulating states, a fact 
that was already pointed out by Isacsson et al. |9| . How- 
ever, many of the degeneracies appearing in their work 
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FIG. 4: (Color online) Comparison between the Gutzwiller 
approach (solid blue lines) and the Gross-Pitaevskii theory 
(dashed red lines) in a three-dimensional system. The pa- 
rameters were computed for a lattice of depth 15 Er in all 
directions. We fix n/Uoo in the Bose-Hubbard model and 
changed Gto/Uoo- We show comparisons for the condensate 
order parameter (ip x ), onsite angular momenta per particle 
\L\/iit, as well as for the fluctuations AL^/tit- In (a), (d), 
and (g) the strong coupling region was in a Mott state with 1 
atom per site, in (b), (e), and (h) the strong coupling region 
was in a Mott state with 2 atoms per site, and in (c), (f), 
and (i) the strong coupling region was in a Mott state with 3 
atoms per site. Note that we choose a specific Mott insulating 
state with nr = 1 so that it had an angular momentum of 1 
per atom. Since this region is in our approximation strongly 
degenerate, many other choices would have been equally jus- 
tified. 



are actually artifacts of utilizing a harmonic approxima- 
tion. Furthermore, we found that the p-band Mott lobes 
follow roughly the structure for the Mott lobes on the 
lowest band, as depicted in Figs. Q] and [31 This is in con- 
trast with the results of Ref. [9j where the Mott lobes ex- 
tends over larger parameter regimes, and they moreover 
show an anomalous behavior with large variations in the 
sizes of neighboring Mott lobes. This discrepancy seems 
to originate from a factor of 2 missing for the cross terms 
proportional to n x n y , n y n z , and n y n z in their work. 

In Fig. [¥]we compare the Gutzwiller approach and the 
Gross-Pitaevskii approach by showing one component of 
the condensate order parameters, angular momenta per 
particle, as well as the fluctuations of the z-component of 
the angular momentum as a function of 6io/£^oo- I n this 
figure we fixed h/Uqq in the Bose-Hubbard model phase 
diagram and changed 6io/C^oo and for each point com- 
puted the corresponding solution of the Gross-Pitaevskii 
equations with the same density. The fixed values of 
/i/C/oo were chosen in such a way that the starting point 
was in the center of the Mott insulating phase with either 
1, 2, or 3 atoms per site. 

With the exception of fluctuations of the single particle 
per site angular momentum, we can see that outside the 
Mott insulating regions the Gross-Pitaevskii theory can 
quickly predict the value of the condensate order param- 



eters quite accurately. Angular momenta are in a sense 
sometimes even better predicted by the Gross-Pitaevskii 
theory, since in the Mott phase with 2 atoms per site an- 
gular momentum is non-zero and behaves qualitatively in 
the same way in the two different approaches. However, 
with 3 atoms per site the angular momentum vanishes in 
the Gutzwiller approach, but is non-zero in the Gross- 
Pitaevskii approach. Also the fluctuations of angular 
momentum agree well in the SF regime. These results 
give us a benchmark for the reliable use of the Gross- 
Pitaevskii formalism for the description of the excited 
band bosons, and we especially find that the mean-field 
treatment is surprisingly accurate even relatively close 
to the Mott boundaries were quantum fluctuations are 
known to become significant. 

Earlier we pointed out the possibility of the broken per- 
mutational symmetry when g xx > 3g xy . In this case the 
order parameters can be unequal. Interestingly, we find 
that this broken symmetry is also reflected in the Mott 
insulating state, where the exact ground state (with two 
atoms per site in this example) carrying angular momen- 
tum changes into a superposition 



fee*|200) + V ^e^|020) 



/p z e 



i<t>3 



002) 



with possibly unequal number of atoms in different fla- 
vors, in contrast to the symmetric state (I28|) . In par- 
ticular, for the state (|29|) the angular momentum van- 
ishes. As one moves to the superfluid phase from the 
Mott phase, the permutational symmetry breaking can 
manifest itself by a single non-vanishing order parame- 
ter (ijj a ) followed by a transition into a state with two 
non- vanishing (and equal) order parameters @. 



V. CONCLUSIONS 

In this paper we have explored the properties of 
bosonic atoms on the first excited band of an optical 
lattice. By computing the phase diagrams for two- and 
three-dimensional systems, we found Mott-insulating and 
superfluid phases with more subtle quantum properties 
than those appearing in the lowest band Hubbard model. 
Furthermore, we compared the Gutzwiller theory to the 
Gross-Pitaevskii approach and established the parame- 
ter regimes where the latter description provides a good 
approximation to the physical system. 

Here we found that bosons on the p-band can 
form a staggered-vortex superfluid composed of anti- 
ferromagnetically ordered vortices and anti-vortices. Ro- 
tation breaks the degeneracy of the vortex and anti- 
vortex state and it would be interesting to explore how 
rotation favoring vortex lattice formation competes with 
the physics of staggered- vortex superfiuids. Also, in fairly 
shallow lattices where effects due to interactions can be 
pronounced, dispersions can develop swallowtails in the 
vicinity of the Brillouin zone edge and period doubled 
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states can appear [29|, [3(| ■ In the previous analysis which 
assumed an one-dimensional systems, the swallowtails 
were found to be related to the existence of solutions 
corresponding to a train of solitons. It would be of inter- 
est to explore the similar situation in higher dimensions, 
where stability properties are often very different. 

Experiments are typically done in optical lattices with 
an additional trapping potential acting at the back- 
ground. Here we studied only the homogeneous solutions 
and this assumption is valid locally when the background 



trapping potential varies slowly compared to the lattice 
spacing. Our results can be applied in a trap using the 
local density approximation or by adding a site depen- 
dent energy offset to the Hamiltonian. However, the size 
of computations using the multi-flavor Gutzwiller ansatz 
grow quickly as a function of system size which at this 
stage limits us to fairly small systems. Inhomogeneous 
density distribution is easier to take into account within 
the mean-field approximation. 
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